#ifndef UTILS_CC
#define UTILS_CC

#include<string>
#include<fstream>
#include<sstream>
#include<iostream>
#include<algorithm>
#include<iterator>
#include<vector>
#include<map>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include "def.hpp"
#include "utils.hpp"
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_randist.h>


using namespace std;

double get_sum(double*& buffer, int& K){
  double result = 0.0;
  for(int k=0; k < K; k++){
    result += buffer[k];
  }
  return result;
}

double log_sum(double log_a, double log_b){
  if(log_a < log_b)
    return log_b + log(1 + exp(log_a - log_b));
  else
    return log_a + log(1 + exp(log_b - log_a));
}

double log_sub(double log_a, double log_b){
  if(log_a < log_b)
    return log_b + log(1 - exp(log_a - log_b));
  else
    return log_a + log(1 - exp(log_b - log_a));
}

double log_sum_array(double*& buffer, int K){
  // first find the maximum
  double max_x = buffer[0];
  double result = 0.0;  
  for(int k=0; k < K; k++){
    if(max_x < buffer[k])
      max_x = buffer[k];
  }
  // update
  for(int k=0; k < K; k++){
    result = result + exp(buffer[k] - max_x);
  }
  result = log(result) + max_x;
  return result;
}

double log_sum_vector(vector<double>& buffer){
  // first find the maximum
  double max_x = buffer[0];
  double result = 0.0;
  for(unsigned int k=0; k < buffer.size(); k++){
    if(max_x < buffer[k])
      max_x = buffer[k];
  }
  // update
  for(unsigned int k=0; k < buffer.size(); k++){
    result = result + exp(buffer[k] - max_x);
  }
  result = log(result) + max_x;
  return result;
}


double digamma(double x){
  return gsl_sf_psi(x);
}

double log_gamma(double x){
  return gsl_sf_lngamma(x);
}

double sum_log_gamma(double*& buffer, int& K){
  double s = 0.0;
  for(int i=0; i < K; i++){
    s = s + log_gamma(buffer[i]);
  }
  return s;
}

double random_gamma_distribution(double shape, double scalar){
  return gsl_ran_gamma(PuzaDEF::Instance()->gBaseRand, shape, scalar);
}

double random_uniform_distribution(double lower, double upper){
  return gsl_ran_flat (PuzaDEF::Instance()->gBaseRand, lower, upper);
}

double l1_distance(double*& a, double*& b, int& K){
  double result = 0.0;
  for(int k=0; k < K; k ++){
    result = result + fabs(a[k] - b[k]);
  }
  return result;
}

double l2_distance(double*& a, double*& b, int& K){
  double result = 0.0;
  for(int k=0; k < K ; k++){
    result = result + (a[k] - b[k]) * (a[k] - b[k]);
  }
  result = sqrt(result);
  return result;
}

double beta_function(double a, double b){
  return gsl_sf_beta(a, b);
}

double beta_distribution(double a, double b, double x){
  return gsl_ran_beta_pdf(x, a, b);
}

double Min(double a, double b){
  if(a < b)
    return a;
  else
    return b;
}

double Max(double a, double b){
  if(a > b)
    return a;
  else
    return b;
}

double Trunc(double vj, double alpha){
  if(vj <=0)
    return Min(0,vj + alpha);
  else
    return Max(0,vj - alpha);
}


#endif
